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ABSTRACT 
Control  systems  containing  two  nonlinear  elements  are  separated 
into  three  broad  classes  which  are  dependent  upon  the  location  of  an 
initial  signal  excitation  with  respect  to  the  two  nonlinear  elements. 
The  nonlinear  elements  are  considered  as  the  "<■  and  Q  variables  of  the 
parameter  plane  which  defines  the  region  containing  the  roots  of  the 
system  characteristic  equation.  A  "dynamic  describing  function" 
curve,  which  determines  the  locus  of  the  M  or  operating  point  for  in- 
itial condition  input  signals,  is  plotted  on  the  parameter  plane. 
Additionally,  the  transient  response  is  graphically  determined  for 
conditions  of  second  order  dominance.   The  theoretical  results  are 
compared  with  actual  analog  computer  simulation. 
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CHAPTER  I 
INTRODUCTION 

The  significant  feature  of  nonlinear  control  problems  is  that  no 
two  problems  are  identical  --  each  one  has  a  unique  solution.   At  best 
one  can  discuss  classes  of  nonlinear  problems  where  the  nonlinear  ele- 
ments appear  in  the  error  channel,  the  feedback  loop,  or  in  combinations 
in  the  error  channel  and  the  feedback  channel.   In  keeping  with  this 
qua3i-clas8  criterion,  control  systems  containing  two  amplitude -sensi- 
tive nonlinearities  which  are  related  with  respect  to  a  common  signal 
excitation  will  be  covered  in  this  thesis. 

Control  systems  containing  single  nonlinear  elements  have  been  ac- 
curately analyzed.   Siljak  [5]  has  used  the  parameter  plane  technique 
for  the  analysis  and  synthesis  of  single  nonlinear  problems  with  ex- 
cellent results.   In  Wilson's  thesis  [8]  the  limit  cycles  of  feedback 
control  problems  were  determined  by  Mitrovic's  method  and  checked  by 
the  root  locus  plot,  Nyquist  diagram,  and  analog  simulation  with  high 
correlation.   Another  technique  which  has  been  used  to  predict  the 
transient  performance  to  within  five  to  ten  percent  accuracy  for  a  step 
input  signal  is  MacAllister's  method  [2]  of  calculating  a  "gain  phase" 
function  for  the  nonlinearity  based  on  a  damped  cosinusoidal  input 
signal.   An  obvious  brute  force  method  of  solving  the  nonlinear  problem 
is  to  use  linear  piecewise  approximations  of  the  nonlinearity  and  pro- 
gram the  digital  computer  to  handle  the  iterative  calculations.  Al- 
most any  degree  of  accuracy  can  be  obtained  to  get  the  transient  response 
by  digital  computer  solution  of  the  system's  differential  equation. 


For  systems  containing  two  nonlinearities  the  determination  of 
system  performance  has  been  somewhat  less  successful  due  to  the  in- 
creased complexity  of  the  problem.   One  of  the  most  significant  contri- 
butions is  Siljak's  paper  [5]  on  the  use  of  the  parameter  plane  for 
establishing  the  stability  region  and  limit  cycles  of  nonlinear  control 
systems  with  respect  to  the  system's  parameters  and  initial  conditions. 
The  parameter  plane  technique  will  be  used  in  this  thesis  for  determin- 
ing stability  and  limit  cycles  and,  in  addition,  as  a  tool  for  calculat- 
ing the  transient  response  under  conditions  of  second  order  dominance. 
An  unpublished  paper  by  Yockey  and  Thaler  [9]  obtains  the  transient 
response  of  a  nonlinear  system  containing  two  saturation  elements  in  a 
multiloop  feedback  path.   Excellent  correlation  with  analog  computer 
simulation  was  obtained  in  determining  the  transient  response  of  their 
system.   The  concept  of  a  "dynamic  describing"  function  introduced  in 
this  last  paper  will  also  be  used  in  this  thesis. 

While  many  nonlinear  problems  can  be  solved  using  a  digital  com- 
puter or  by  simulation  on  an  analog  computer,  a  control  designer  would 
like  to  know  a  priori  what  will  happen  to  a  system  when  certain  para- 
meters in  a  nonlinear  control  problem  are  varied.   One  design  philosophy 
currently  used  in  control  engineering  is  the  perturbation  technique  where 
small  changes  in  system  parameters  are  made  and  the  system  response  cal- 
culated on  the  digital  computer.   After  several  and  sometimes  many  alter- 
ations to  the  system  parameters  the  desired  response  is  obtained.  While 
this  method  is  obviously  very  accurate  in  the  final  analysis,  consider- 
able time  is  wasted  in  the  design  phase.   Furthermore,  turning  the  dials 
of  an  analog  computer  or  submitting  new  data  to  a  digital  computer  to 


obtain  a  desired  response  will  not  answer  the  fundamental  questions 
of  "why"  or  "how"  will  the  system  perform. 

In  this  thesis  precise  answers  for  nonlinear  control  systems  will 
not  be  sought.   Rather,  the  stability  regions,  the  limit  cycle  frequency, 
and  the  form  of  the  transient  response  will  be  predicted  for  two  classes 
of  nonlinear  systems  and  discussed  for  a  third  class. 

On  the  basis  of  the  work  done  by  Siljak  [5]  it  appears  that  two 
nonlinear  elements  can  be  conveniently  classified  with  respect  to  an 
initial  signal  excitation  into  three  general  cases. 

Case_l.   Identical  signal  excitation  to  both  nonlinear  elements. 

In  Figure  la  the  signal  denoted  by  X  is  the  same  input  signal  to 
nonlinear  elements  N  and  N  .   The  transfer  functions  denoted  by  G(s) 
are  assumed  to  provide  sufficient  low  frequency  filtering  for  the  con- 
trol system.   This  example  is  considered  the  simplest  case. 

Case  2.   Input  signal  excitation  to  nonlinear  elements  related  by  a 
linear  differential  equation. 

This  case  is  shown  in  Figure  lb  where  the  input  signals  to  N.  and 
N„  are  respectively  X  G  . (s)  and  X.   The  input  signal  to  N  is  a  function 
of  frequency  and  amplitude  while  the  input  signal  to  N  is  a  function  of 
amplitude  only.   The  additional  frequency  sensitivity  of  this  system 
complicates  the  problem. 

Case  3.   Input  signal  excitation  to  nonlinear  elements  related  by  a 
nonlinear  differential  equation. 

The  most  difficult  case  appears  in  Figure  lc  where  the  signal  in- 
put to  N  is  X  and  the  input  signal  to  N  is  x9  N  1  £g  (S))  where  the 
brackets  are  intended  to  represent  some  transform  relationship  rather 
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Figure  1.   General  Classification  of  Control  Systems 
with  Two  Nonlinear  Elements, 


than  interpreting  the  three  quantities  in  a  strict  multiplicative  sense, 

These  three  cases  will  be  covered,  the  first  two  experimentally 
and  the  last  case  theoretically. 


CHAPTER  II 
APPLICATION  OF  THE  PARAMETER  PLANE  TO  THE 
TWO  NONLINEARITIES  PROBLEM 
The  method  used  to  determine  system  performance  in  this  thesis  is 
based  on  the  superposition  of  a  dynamic  describing  function  curve  on 
the  parameter  plane  curves  for  a  control  system  containing  two  ampli- 
tude sensitive  nonlinearities.   This  chapter  reviews  the  theory  behind 
this  analytical  technique. 

2-1  The  parameter  plane  technique 

In  the  previously  mentioned  reference  by  Siljak  it  has  been  shown 
that  the  sum  of  the  terms  of  the  characteristic  equation  can  be  represent- 
ed as  a  sum,  i.e. , 

f(s)  =   £  a  sk  (2-1) 

K.-0      K 

where  the  a,  are  the  linear  coefficients  of  the  characteristic  equation. 


By  letting  s  =  -£,uJr.+   j  ,Av\|  I  "  I  equation  (2-1)  can  be  expanded 

in  terms  of  Chebyshev  polynomials  of  the  second  kind  which  are  func- 

p 
tions  of  6  ,  the  damping  ratio. 

The  coefficients  of  the  characteristic  equation  can  be  expressed 

as 

■k"  t*-\  +  S<k  +  \  <2-2> 

where  alpha  and  beta  are  any  variable  system  parameters.   This  feature 
of  the  parameter  plane  technique  gives  the  analyst  or  designer  of  con- 
trol systems  two  degrees  of  freedom  as  opposed  to  the  single  parameter 
techniques  such  as  root  locus,  Bode  diagram,  or  Nyquist  plot.   Since  an 
amplitude  sensitive  nonlinearity  can  be  considered  as  a  variable  gain 


element  whose  value  varies  between  zero  and  one,  it  is  convenient  to 
let  alpha  and  beta  become  the  variable  nonlinear  elements  when  the 
characteristic  coefficients  are  expressible  in  the  form  of  equation  (2-2). 

Although  the  two  parameters  alpha  and  beta  can  be  solved  for  alge- 
braically from  equation  (2-1),  it  is  more  convenient  to  plot  a  net  of 
curves  in  terms  of  £   and  (xj^  ,  J>  and  10  ,  or  i^r,,  and  uJ  ,  depending 
upon  what  particular  system  characteristics  are  under  consideration. 
In  the  following  discussion  the  utility  of  the  different  nets  of  curves 
will  be  discussed. 

2-1-1     Zeta  -  omega  curves 

Because  of  the  recursive  nature  of  the  Chebyshev  polynomials  it  is 
convenient  to  plot  a  net  of  curves  for  the  variables  alpha  and  beta  in 
terms  of  J>  ,  the  damping  ratio,  and  U-^  ,  the  undamped  natural  frequency. 
A  computer  program  which  will  plot  this  net  of  curves  appears  in  Nutting's 
thesis  [3].   Since  engineers  are  used  to  the  concept  of  the  numeric,  /  , 
the  utility  of  such  a  set  of  curves  is  obvious.   The  controls  engineer 
can  tell  at  a  glance  what  effect  a  change  in  one  or  both  of  the  system 
variables  will  have  on  system  stability  by  the  location  of  the  M  or 
operating  point  with  respect  to  the  J   =  0  curve.   In  addition  to  deter- 
mining the  relative  stability,  the  presence  or  absence  of  possible  limit 
cycles  can  be  established.   Additionally,  if  second  order  dominance  is 
justified  an  indication  of  the  transient  response  can  be  determined  since 
the  complex  roots  are  obtained  directly  from  the  graph  by  reading  the 
values  of  J   and  0JK  at  the  operating  point.   Any  real  roots  may  be  ap- 
proximated by  drawing  tangents  to  the  J    =  1  curve. 

The  previous  discussion  considered  a  static  case  or  a  linear  system. 


For  nonlinear  systems  where  alpha  and  beta  are  the  two  nonlinear  ele- 
ments, a  dynamic  operating  condition  exists  since  the  nonlinearities 
are  not  only  sensitive  to  the  input  signal  amplitude  but  also  to  fre- 
quency.  To  cover  this  nonlinear  situation,  describing  function  theory 
will  be  discussed  in  a  later  section. 

Since  describing  function  theory  depends  upon  the  oscillating  fre- 
quency, uj  ,  it  is  necessary  to  modify  the  constant  *M\  curves  to  con- 


stant UJ  curves  by  the  relationship  ^  =  ^M-\N  '  ~  I         .      The  recursive 
relationships  developed  using  Chebyshev  polynomials  are  not  amenable  to 
this  direct  substitution.   However,  a  set  of  constant  ^0  curves  can  be 
obtained  from  the  binomial  coefficient  expansion  described  in  the  next 
section.   Hence,  it  is  possible  to  obtain  a  net  of  J~0j     curves  which 
may  be  used  directly  with  describing  function  theory. 

2-1-2     Sigma  -  omega  curves 

The  concept  of  constant  zeta  curves  or  constant  omega  curves  is 
accepted  in  engineering  circles.   But  the  utility  of  constant  sigma  or 
t,^  curves  is  not  generally  recognized  since  there  has  been  no  real 
need  for  them.   Perhaps,  the  most  difficult  idea  behind  constant  sigma 
curves  is  that  engineers  have  no  intuitive  feeling  for  this  variable. 
However,  the  reference  curve  for  the  sigma  curves  is  the  sigma  =  0  curve 
which  is  identical  to  the  zeta  =  0  curve.   The  usefulness  of  the  con- 
stant sigma  curves  will  become  apparent  when  the  transient  response  of  a 
system  is  desired.   In  those  cases  where  second  order  dominance  is  justi- 
fied, which  is  usually  the  case  for  well-designed  control  systems,  the 
value  of  sigma  determines  the  rate  of  exponential  signal  decay  in  the 
stable  region. 
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Recently,  Siljak  has  shown  that  the  individual  terms  comprising  the 
characteristic  equation.,  i.e. ,  f(s)  =    ^   as   can  be  expanded  in 
terms  of  binomial  coefficients.   Consider  the  rectangular  coordinate  re- 


lationship 


=k  ■  \ +  j\  <2"3> 


where  X  =  X  (  ^  oj  )  and  Y  =  Y  (  0",  uj  )  .   The  first  order  term  of 
equation  (2-3)  can  be  expressed  as  s  =  X  +  jY  ,  or  s  =  C+  j  uJ  .   If 
the  latter  equation  is  substituted  into  equation  (2-3)  the  following 
pair  of  recursion  formulas  are  obtained  by  induction: 

*k+l    "   2Xl\  +   <V  +  Yl2)     Xk-1   =  °  (2-4) 

\+i  -  "A  +  <xi2  +  Yi2)  \-i =  ° 

subject  to  the  following  initial  conditions 


X0  =  l  Y0  "  ° 


x  -  <r  y  -  oJ 


(2-5) 


Each  of  the  k+1  terms  in  equations  (2-4)  is  a  binomial  coefficient 
which  satisfies  the  following  pair  of  equations 

1  aA  "  °  (2-6) 

to  Vk  =  ° 

The  first  equation  in  (2-6)  represents  the  summation  of  the  real  terms 
and  the  second  equation  the  summation  of  imaginary  terms.   Alpha  and 
beta  of  the  linear  coefficient  a  =  b,<*-  +  c  (3  +   d   are  now  implicit 
functions  of  sigma  and  omega?i.e.,  <*"-  =  <X.  (  tf)   uJ  )  and   ^  =   @  (  (J^  uj)  . 
The  following  pair  of  transform  equations  developed  in  Siljak' s  original 
paper  [5]  still  holds: 

Bl^+Cl^+Dl-°  (2-7) 

B2=c  +  C2  @    +  D2  -  0 


where  the  capitalized  variables  are  implicit  functions  of  sigma  and  omega 
and  are  equal  to  the  summation  of  the  corresponding  lower  case  lettered 
terms  of  the  generalized  coefficient  a,  =  b  °<-  +  c  (3     +  d  .   From 
Cramer's  rule  alpha  and  beta  can  be  solved  for  from  equations  2-7: 

co  -   C1D2  -  C2D1 

B1C2  "  B2C1 
1  l  l  l  (2-8) 

Q     -       B2D1  -  B1D2 
B1C2  "  B2C1 

Equations  (2-8)  are  readily  programmed  on  the  digital  computer  with  the 

aid  of  recursion  formulas  (2-4).   A  computer  program  which  plots  a  net 

of  sigma  -  omega  curves  for  the  parameters  alpha  and  beta  appears  in 

Appendix  B.   As  long  as  the  coordinate  axes  dimensions  for  alpha  and 

beta  remain  constant  the  constant  omega  curves  of  this  section  can  be 

superimposed  on  the  constant  zeta  curves  of  the  previous  section  to 

obtain  a  net  of  constant  S~(J^     curves.   These  two  nets  of  curves, 

namely,  '[  -  WJ  and   0~-CO  ,  will  be  used  in  subsequent  discussions. 

Since  1      is  the  more  generally  understood  term,  it  will  be  used  more 

frequently  than  sigma. 

2-2  Describing  function  technique 

In  the  analysis  of  nonlinear  systems  it  is  necessary  to  make  certain 
linear  approximations  concerning  the  nonlinear  elements  to  obtain  some 
practical  idea  of  system  performance.   One  of  these  linearization  methods 
is  the  describing  function  in  which  a  sinusoidal  signal  of  the  form  A 
sin  CO  t  is  applied  to  the  nonlinearity .   If  the  input  -  output  character- 
istic of  the  nonlinear  element  is  known,  the  output  can  be  expanded  in  a 
Fourier  series.   The  ratio  of  the  fundamental  component  of  the  Fourier 

10 


expansion  to  the  peak  amplitude  of  the  sinusoidal  input  signal  is  de- 
fined as  the  describing  function,  or  more  explicitly  the  gain-phase 
characteristic  of  the  nonlinear  element. 

2-2-1     Inherent  shortcomings 

The  assumption  of  a  single  sinusoidal  input  signal  of  constant 
amplitude  occurs  in  a  nonlinear  system  only  under  oscillatory  condi- 
tions,^.e.,  J     =  0.   Furthermore,  the  assumption  of  a  single  frequency 
component  in  the  describing  function  derivation  is  valid  only  if  suf- 
ficient low  frequency  filtering  is  available  in  the  control  system  to 
limit  the  harmonics  of  the  limit  cycle  to  negligible  values.   For  a  sec- 
ond order  system  this  assumption  would  normally  not  be  valid,  but  for 
higher  order  systems  the  accuracy  of  the  describing  function  technique 
would  normally  improve.   In  practice  the  describing  function  gives  the 
best  results  for  third  order  systems  or  higher  under  limit  cycle  con- 
ditions. 

Since  the  nonlinear  element  is  a  function  of  the  input  signal,  A, 
and  the  oscillating  frequency,  OJ  ,  the  only  region  on  the  parameter 
plane  where  the  describing  function  is  justified  occurs  when  UJ  =  ^yn 
on  the  «i  =  0  or   0=0  curve.   But  the  operating  point,  M(  °<-  ,  Q  ) 
or  M(N  ,N  )  where  the  alpha  -  beta  parameters  become  the  nonlinear  vari- 
ables N  and  N  ,  moves  in  the  parameter  plane  due  to  the  dependence  of 
the  nonlinearity  on  A,  .[  ,  and  Cu^  .   Consequently,  it  is  necessary  to 
modify  the  constant   U-)^  curves  to  constant  uJ  curves  by  the  relation- 
ship °J     -    °^r^{  I  "  i      .   In  the  previous  section  a  recursion  formula 
using  binomial  coefficients  was  indicated  to  obtain  constant  tO  curves 
which  can  be  used  with  the  variables   vi  or  (f   . 
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For  two  nonlinearities  where  one  or  both  are  a  function  of  fre- 
quency, the  loci  of  the  describing  functions  will  be  constant  <0  lines 
for  varying  signal  input.   The  intersections  of  the  common  radian  fre- 
quencies between  the  parameter  plane  curves  and  the  describing  func- 
tion curves  will  determine  the  values  of  the  alpha  -  beta  variables 
for  a  particular  initial  condition  input.   A  series  of  these  points 
connected  by  a  smooth  curve  is  called  a  "dynamic  describing"  function 
curve  in  the  aforementioned  paper  by  Yockey  and  Thaler.   In  summary 
this  means  that  for  a  specified  initial  condition  the  operating  or  M( 
*us   ,  §   )  point  should  lie  on  the  dynamic  describing  function  locus. 
The  superposition  of  the  dynamic  describing  function  locus  atop 
the  parameter  curves  places  in  evidence  the  roots  of  the  characteris- 
tic equation  at  each  instant  of  time.  The  motion  of  the  M(  <^  ,  §  ) 
point  is  due  to  the  changed  values  of  the  nonlinear  elements.   This 
important  information  can  be  used  to  obtain  the  transient  response  of 
the  system  under  second  order  dominancy  conditions,  i.e.,  where  one 
pair  of  complex  poles  is  dominant.   A  brief  discussion  of  the  theory 
and  an  example  using  J  -  oJ   curves  and   <T  -  ^  curves  is  worked  in 
Appendix  A. 

2-2-2     Application  of  the  "dynamic  describing  function"  locus 

In  view  of  the  previous  discussion  a  mathematical  justification  for 
the  describing  function  locus  for  two  nonlinearities,  which  are  related 
in  one  of  the  three  ways  indicated  in  the  introduction,  will  not  be  at- 
temped.   Rather,  the  problem  will  be  looked  at  from  an  engineering 
viewpoint.   Other  system  variables  such  as  phase  shift,  real  root  ef- 
fect, and  complex  roots  of  higher  order  systems  may  alter  physical 
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reasoning  as  the  dynamic  describing  function  locus  gets  further  away 
from  either  the  I  ■  0  or  the  <T  =  0  curves.   However,  these  effects 
are  not  apparent  intuitively,  let  alone  mathematically.  Using  this 
semi-empirical  approach,  the  amplitude  of  the  reference  signal  excita- 
tion, A,  will  be  considered  as  the  peak  amplitude  of  a  sinusoidal 
signal  at  one  instant  of  time. 
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CHAPTER  III 
TWO  NONLINEAR  ELEMENTS  WITH  IDENTICAL  SIGNAL  EXCITATION 

The  first  case  of  two  nonlinear  elements  related  by  a  common  sig- 
nal excitation  is  considered  in  this  chapter.   A  third  order  system  con- 
taining two  amplitude  sensitive  elements  is  chosen  as  the  working  model. 
Based  on  the  previous  discussion  in  Chapter  II  on  the  parameter  plane 
technique  and  the  describing  function  method,  the  system  performance  is 
predicted  and  compared  with  analog  simulation. 

3-1  The  working  model 

The  experimental  system  is  shown  in  Figure  3-1.   The  uncompensated 

plant  with  the  nonlinear  elements  is  inherently  unstable,  and  is  stabil- 

zed  by  tachometer  feedback.   With  the  insertion  of  the  nonlinear  elements, 

N-  and  N  ,  the  prediction  of  system  performance  is  complicated.   The  two 

dead  zone  elements  have  a  common  signal  excitation,  A,  which  is  assumed 

to  be  some  initial  condition  voltage.   The  dead  zone  width  denoted  by 

+  E,  does  not  have  to  be  identical  for  both  elements. 
—  d 

The  characteristic  equation  for  this  system  is 

s3  +  10s2  +  (ION  +  ION  )s  +  100  N  =  0  •         (3-1) 

Using  the  parameter  plane  technique  it  is  convenient  to  let   ^  =  N 
and   $  =  N  .   Equation  (3-1)  becomes 

s3  +  10s2  +  (10  <*-  +  10  3  )s  +  100*.=  0  (3-2) 

A  set  of  zeta-omega  and  sigma-omega  curves  can  be  plotted  for  this  sys- 
tem using  digital  computer  programs.   In  addition  the  instantaneous  values 
of  N.  and  N-  can  be  found  for  different  amplitudes  of  input  signals  and 
various  dead  zone  widths.   From  reference  [7]  the  describing  function  for 
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Figure   3-1.      Block  Diagram   of  Third   Order   System 
Mr-?  4-v.    rv.rri   TvTonl -i  nftflT*    Elp.ments. 
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Plgure  3-2.      Dynamic   Describing  Function  Curves   on 
Zeta-Omega   Curves. 
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Plgure  3-3. 


Dynamic   Describing  Function  Curve  on 
Sigma-Omega   Curves. 
Hi   -   Dead   Zone  Ng   -   Dead   Zone 


Ed  -  i  5  Volts 
E<J  :  t  3  Volts 
Ed  =  t  3   Volts 


Ed  -  1  3  Volts 
Ed  =  ±  3  Volts 
Ed  =  ±5   Volta 


dead  zone  is 


Gd  =  IF  (  \    Vin ^  R"R  ^"r2)  (3"3) 


where  R  =  EVA,  E,  being  one-half  '>;he  total  dead  zone  width  and  A  the 
d     d  fj 

peak  amplitude  of  the  signal  input.   This  equation  is  readily  programmed 
on  the  computer  to  obtain  the  values  of  N-  and  N  for  various  amplitudes 
of  A  and  ratios  of  R.   Plotting  these  parameter  values  on  top  of  either 
the  zeta-omega  or  sigma-omega  curves  determines  the  systems  dynamic  des- 
cribing function. 
3-2  Graphical  results 

Figures  3-2  and  3-3  are  the  zeta-omega  and  sigma-omega  curves  respec- 
tively for  the  system.   Physically,  the  region  of  interest  on  both  sets 
of  curves  is  the  first  quadrant  since  N.  and  N  are  passive  devices  hav- 
ing values  lying  between  zero  and  one.   From  the  zeta-omega  curves,  the 
normal  stable  region  is  further  confined  between  the  I-  0   and  the  J,    =  I 
curve  in  the  first  quadrant.   As  noted  on  both  the  zeta-omega  and  sigma- 
omega  curves,  the  ^or  N  axis  is  the  I  =  0  or  3~  =  0  curve. 

On  top  of  both  the  zeta-omega  and  sigma-omega  curves  are  superimposed 
the  describing  function  loci  for  three  curves:   (1)  dead  zone  width  of 
N  =  +  5  volts  and  the  dead  zone  width  of  N  =  +  3  volts,  (2)  dead  zone 
width  of  both  N  and  N  =  +  3  volts,  and  (3)  dead  width  of  N  =  +  3  volts 
and  dead  zone  width  of  N  =  +  5  volts.   Since  both  nonlinear  elements  are 
functions  of  signal  amplitude  only,  the  smooth  curves  connecting  the 
labelled  points  also  determine  the  dynamic  describing  function  curve  for 
the  three  cases. 

Dynamic  describing  function  curve  1  reaches  the  3=1  curve  before 
the  nonlinear  effects  are  eliminated.   Since  £   is  greater  than  one  when 
the  N?  or  (3  boundary  is  reached,  the  system  response  will  be  overdamped 
as  modified  by  the  nonlinear  element  N9.   Curve  two  should  theoretically 
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Figure   5-4.      Analog   Computer   Simulation, 
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end  in  the  nebulous  dead  zone  region  which  occurs  when  N  =  N  =  0.   In 
the  case  of  curve  3,  the  £   =  0  curve  is  reached  when  the  envelope  of 
the  exciting  signal  reaches  5  volts.   Since  smaller  amplitudes  are  in- 
dicated below  this  point,  a  stable  limit  cycle  is  predicted  for  curve  3. 

3-3  Comparison  with  analog  computer 

The  analog  computer  setup  for  this  two  nonlinear  element  ampli- 
tude sensitive  only  system  appears  in  Figure  3-4.   Initially,  the  com- 
puter was  set  up  for  the  nonlinear  configuration  for  curve  3  to  com- 
pare the  actual  limit  cycle  performance  with  the  dynamic  describing 
function  locus  prediction.   The  results  are  shown  below. 


Limit  cycle 

Amplitude 

uj 

Theoretical 

5.0  volts 

1.67 

Computer 

5.0  volts 

1.65 

To  check  the  system  transient  response  an  initial  condition  of  11 
volts  was  assumed  for  each  of  the  dynamic  describing  function  curves  in 
Figures  3-2  and  3-3.   The  transient  response  was  calculated  using  the 
sigma-omega  curves  as  demonstrated  in  Appendix  A.   These  results  ap- 
pear in  Figures  3-5  to  3-7.   In  the  upper  right  hand  corner  is  the 
actual  brush  recording  trace  of  the  system  response  for  the  indicated 
initial  condition  and  nonlinear  configuration. 

The  transient  response  appearing  in  Figure  3-5  for  curve  1  reason- 
ably approximates  the  analog  results  until  the  signal  envelope  reaches 
5  volts.   Theoretically,  the  response  should  be  some  slowly  decaying 
waveform  because  of  the  high  values  of  X   .   Practically,  the  system 
had  an  extra  overshoot  and  then  a  long  decaying  response.   For  this 
unique  case  part  of  the  shortcoming  lies  in  the  graphical  technique 
which  is  limited  to  damped  oscillatory  systems.   For  curve  two  the 
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transient  prediction  in  Figure  3-6  is  improved  because  the  high  values 
of  zeta  are  no  longer  prominent.   During  the  first  five  seconds  of  the 
transient  response  the  theoretical  results  closely  approximate  the  ana- 
log results.   By  this  time  the  signal  envelope  is  less  than  three  volts 
and  the  system  is  in  the  nebulous  dead  zone  region  where  NL  ■  N_  =  0. 
The  analog  computer  indicates  some  non-oscillatory  waveform  after  this 
terminal  point  is  reached.   The  best  correspondence  between  the  theoreti- 
cal calculations  and  the  analog  simulation  is  shown  in  the  transient  re- 
sponse curves  appearing  in  Figure  3-7.  Two  overshoots  are  indicated  on 
both  the  calculated  transient  response  and  the  actual  analog  response  be- 
fore the  system  reaches  the  stable  limit  cycle.   As  indicated  earlier  the 
amplitude  and  frequency  of  the  limit  cycle  agreed  very  closely. 

3-4  Summary  of  results 

This  example  indicates  that  the  effect  on  the  variations  in  system 
performance  due  to  changes  in  the  nonlinear  variables  can  be  predicted 
to  some  extent  by  noting  the  changing  dynamic  describing  function  curves. 
As  already  noted  excellent  prediction  of  system  performance  results  when 
the  jf  values  are  low  or  in  the  vicinity  of  the  £,   =  0  curve.   When  the 
dynamic  describing  function  curve  lies  in  the  rapidly  changing  high  zeta 
or  overdamped  regions  of  the  parameter  plane,  it  becomes  increasingly 
difficult  to  anticipate  system  performance  accurately. 
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CHAPTER  IV 

TWO  NONLINEAR  ELEMENTS  RELATED  TO  A  COMMON  SIGNAL 

BY  A  LINEAR  DIFFERENTIAL  EQUATION 

In  order  to  establish  the  system  performance  for  nonlinear  ele- 
ments related  to  a  common  reference  signal  by  a  linear  differential 

equation,  nonlinearities  N  and  N_  for  case  2  are  considered  as  satura- 

1      2 

tion  or  dead  zone  elements.   The  theoretical  describing  function  curve 
is  calculated  for  several  test  cases.   A  graphical  technique  which 
assumes  second  order  dominance  was  used  to  calculate  the  transient  re- 
sponses.  The  systems  were  simulated  on  the  analog  computer  and  compari 
sons  were  made  between  the  actual  and  theoretical  amplitudes  and  fre- 
quencies of  the  limit  cycles  as  well  as  the  system  transient  responses. 

4-1   Saturation  in  the  error  channel  and  dead  zone  in  the  tachometer 

feedback  channel 

Figure  4-1  shows  the  working  model  with  the  dead  zone  width  of  N 
and  the  saturation  level  of  N  set  at  +  1  volt.   Using  the  graphical 
method  outlined  in  section  2-2-1  a  dynamic  describing  function  curve, 
which  appears  in  Figures  4-2  and  4-3,  was  determined  for  this  system. 
As  indicated  in  section  2-2-1  the  loci  of  the  describing  function 
curves  are  superimposed  on  the  zeta-omega  and  sigma-omega  parameter 
plane  curves.   The  small  circles  indicate  the  points  of  common  radian 
frequency,  (A.)  ,  between  the  describing  function  loci  and  the  parameter 
plane  curves.  Adjacent  to  these  solution  points  are  indicated  the 
amplitude  theoretically  required  for  initial  condition,  A.   The  common 
radian  frequency  is  the  constant  omega  curve  passing  through  the  in- 
dicated points.   The  dynamic  describing  function  curve  is  denoted  by 
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Figure  4-2.   Dynamic  Describing  Function  Curve  on 
Zeta-Omega  Curves. 
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Figure  4-3.      Dynamic   Describing  Function  Curve  on 
Sigma-Omega   Curves. 
Hi  -  Dead  Zone  B^  =  ± 1  Volt 

Kg  -   Saturation  BMt  =  t  1  Volt 


Figure  4-4.  oj  -A  Orld  for  Figures   4-2   and  4-3. 
«!   -   Dead  Zone  Ed  =  ±   1  Volt 

No   -   Saturation  E«at  =  ±     1  v°l* 


the  series  of  heavy  dots.   Figure  4-4  is  the  digital  computer  generated 

(O  -A  grid  for  the  describing  function  loci.   If  tracing  paper  were 
used,  this  grid  could  be  placed  atop  the  modified  parameter  plane 
curves,  and  the  common  radian  frequency,  (O  ,  and  initial  condition 
amplitude,  A,  obtained  by  inspection.   Since  the  parameter  plane  curves 
are  for  any  two  amplitude  sensitive  non-linearities,  it  is  only  neces- 
sary to  obtain  a  new  grid  for  a  different  choice  of  nonlinear  elements. 

The  dynamic  describing  curve  in  Figures  4-2  and  4-3  moves  progres- 
sively leftward,  towards  decreasing  values  of  N-  and  increasing  values 
of  N-,  while  the  required  amplitude  of  the  initial  condition,  A,  de- 
creases until  the  linear  M  -  point  is  reached  which  occurs  at  N.  =0 
and  N  =  1.  At  this  point  the  outer  feedback  loop  will  have  no  effect 
on  system  response  while  the  inner  feedback  loop  acts  as  a  linear  system. 
Also  indicated  on  the  dynamic  describing  function  curve  is  the  amplitude 
and  frequency  of  the  limit  cycle  which  occurs  on  the  J.  =  0  and  (T  =  0 
curves. 

The  system  was  simulated  on  the  analog  computer  and  the  output 
response  recorded  for  various  amplitudes  of  A.   Figure  4-5  shows  the 
analog  computer  setup.   The  amplitude  of  the  limit  cycle  was  determined 
by  averaging  the  values  of  the  first  oscillation  for  initial  condi- 
tions of  13.0  and  13.2  volts.   This  problem  indicated  that  a  finite 
amount  of  time  was  required  for  the  analog  computer  solution  to  stabi- 
lize due  to  the  energy  storage  requirements  of  the  transfer  functions. 
The  amplitude  and  radian  frequency  of  the  limit  cycle  are  compared  be- 
low between  the  graphical  method  and  the  computer  results.   The  oscil- 
lating frequency  is  predicted  accurately,  but  the  theoretical  amplitude 
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is  0.48  volts  low 


Limit  Cycle 

Amplitude 

u^> 

Theoretical 

11.32  volts 

1.18 

Computer 

11.8  volts 

1.18 

To  obtain  a  further  comparison  of  the  reliability  of  the  dynamic 
describing  function  method  the  transient  response  was  graphically  deter- 
mined.   As  previously  noted  this  graphical  method,  developed  by  Siljak, 
assumes  second  order  dominance,  meaning  that  the  results  are  only  valid 
for  a  system  with  a  dominant  pair  of  complex  poles.   Hence,  at  best  the 
results  may  have  only  limited  validity  for  a  fourth  order  system  over 
the  dynamic  describing  function  curve.   An  example  of  this  transient 
response  technique  using  both  the  zeta-omega  and  sigma-omega  curves 
appears  in  Appendix  A. 

Figures  4-6,  4-8,  and  4-10  are  the  graphically  calculated  tran- 
sient responses  for  A  =  9,  6,  and  3  volts  respectively.   The  odd-numberp 
ed  figures  4-7,  4-9,  and  4-11  are  the  corresponding  analog  computer 
transient  responses.   In  the  three  computer  cases  it  should  be  noted 
that  a  finite  amount  of  time  is  required  for  the  exponential  decay  to 
begin  after  the  switch  is  turned  on  whereas  the  graphical  solution 
shows  an  instantaneous  decay.   This  finite  lag  time  due  to  imperfect 
simulation  on  the  analog  computer  causes  the  theoretical  results  to  have 
earlier  zero  crossings  and  larger  overshoot  amplitudes.   In  addition  a 
real-root  effect  is  noticeable  on  the  analog  results  shown  by  the  slow- 
ly decreasing  steady  state  output  which  is  especially  apparent  for  A  = 
3  volts.   As  noted  earlier  the  graphical  technique  makes  no  provisions 
for  real-root  effects. 
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Figure  4-12.   Dynamic  Describing  Function  Curve  or 
Zeta-Omega  Curves. 

N}.  -  Dead  Zone      E<j  =  t  3  Volts 
Ng  -  Saturation     E3at  =  ±  1  Volt 
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Figure  4-13.      Dynamic   Describing  Punotlon  Curve  on 
Sigoa-Ctoega  Curvea. 
Ki   -   Dead  Zone  Ed  =  ±   5  Volts 

»2  -   Saturation  B8at  =  ±    1  v°l* 
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Figure  4-14.  oj  -A  Orid  for  Figures  4-12  and  4-13. 
Hi  -   Dead   Zone  E<j   =  ±    S  Volts 

H2   -   Saturation  E3at  =  ±    1   Volt 


From  this  example  one  would  conclude  that  for  this  system  to  ex- 
pect an  exponentially  damped  response  for  initial  condition  signals  of 
less  than  about  11  volts.   For  A  greater  than  11  volts  both  the  dynamic 
describing  function  curve  and  the  computer  indicate  unstable  conditions. 

4-1-1     Effect  of  changing  one  nonlinearity 

When  the  dead  zone  width  was  increased  to  +  3  volts,  the  dynamic 
describing  function  curye  changed.   Figures  4-12  and  4-13  show  the  modi- 
fied dynamic  describing  function  curve  superimposed  on  the  zeta-omega 
and  sigma-omega  parameter  plane  curves  respectively.   The  oj  -A  grid 
for  this  system  appears  in  Figure  4-14.   This  latter  plot  is  the  only 
additional  labor  required  since  the  parameter  plane  curves  do  not  change 
for  the  same  two  nonlinearities  configuration. 

As  in  the  previous  case  the  amplitude  and  frequency  of  the  limit 
cycle  were  compared  between  the  dynamic  describing  function  curve  and 
the  analog  computer. 


Limit  Cycle 

Amplitude 

UJ 

Theoretical 

21.4  volts 

1.13 

Computer 

25.5  volts 

1.14 

In  this  case  the  disparity  between  the  amplitudes  is  about  twenty 
per  cent.   The  reason  for  this  difference  is  that  describing  function 
theory  is  a  first  order  approximation  for  the  signal  amplitude.,  but  it 
is  a  second  order  approximation  for  the  oscillating  frequency.   Figures 
4-15,  4-17 s  and  4-19  are  computed  transient  responses  for  A  ■  16,  12, 
and  8  volts  respectively.   The  corresponding  transient  responses  for 
the  computer  simulation  appear  in  Figures  4-16,  4-18  and  4-20.   Once 
again  the  real  root  effect  is  apparent  in  the  analog  computer  response. 
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4-2  A  family  of  dynamic  describing  function  curves 

The  previous  two  examples  have  shown  that  system  performance  can- 
not be  predicted  exactly  with  the  tools  at  hand,  but  that  information  is 
available  on  the  dynamic  describing  function  curve  although  the  inter- 
pretation might  be  approximate. 

Since  the  parameter  plane  curves  for  the  working  model  are  for 
any  amplitude  sensitive  nonlinearities,  N  and  N  ,  it  was  decided  to 
reverse  the  nonlinear  elements  -  N  is  now  a  variable  saturation  ele- 
ment and  N   is  a  +  3  volt  dead  zone  element.   Figure  4-21  is  a  block 
diagram  of  this  example. 

By  varying  the  saturation  levels,  it  was  anticipated  that  a  series 
of  dynamic  describing  functions  would  be  obtained  that  would  give  some 
indication  which  areas  of  the  parameter  plane  are  of  the  greatest  util- 
ity in  predicting  the  system  performance  using  the  present  techniques. 
A  series  of  dynamic  describing  function  curves  were  constructed  for 
saturation  levels  of  1,  2,  3,  and  5  volts.   This  composite  plot  ap- 
pears in  Figures  4-22  and  4-23. 

A  few  salient  features  are  readily  apparent  from  these  plots. 
First,  the  amplitude  of  the  initial  conditions,  A,  decreases  in  a  left 
-  to  -  right  manner  except  for  curve  A  where  the  initial  condition  ampli- 
tude appears  to  fluctuate  for  high  values  of  J,     before  settling  down. 
Secondly,  all  the  curves  tend  towards  N  =  1  and  N  =  0  for  smaller 
initial  conditions.   This  result  is  expected  because  for  small  initial 
conditions  the  dead  zone  element  is  effectively  out  of  the  system  while 
the  saturation  element  becomes  linear.   Lastly,  all  the  dynamic  describ- 
ing function  curves  predict  a  stable  limit  cycle.   That  is,  when  the 
curves  cross  the  £   =  0  or  (T  =  0  boundary,  A  becomes  smaller  as  the 
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Figure  4-22.   Family  of  Dynamic  Describing  Function 
Curves  on  Zeta-Omega  Curves. 
N2  Dead  Zone      Ed  =  ±  3  Volts 
A  -  Ni  Saturation     Esat  =  ±  1  Volt 
B  -  Ni  Saturation     E3at  =  ±  2  Volts 
C  -  Ni  Saturation     Esat  =13  Volts 
D  -  Hi  Saturation     Esat  -±  5  Volts 
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Figure  4-23.   Family  of  Dynamlo  Describing  Function 
Curves  on  Slgna-Omega  Curves. 
»2  -  Dead  Zone     Ed  =  i    3  Volts 
A  -  Hi  -  Saturation    E8at=  ±   1  v°lt 
B  -  Hi  -  Saturation    Eeat  =  t  2  Volts 
C  -  Ni  -  Saturation    Eaat  =  t  3  Volts 
D  -  Hi  -  Saturation    Esat  =  ±  6  Volts 
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M-point  enters  the  negative  J     or  unstable  region.   If  the  amplitudes 
were  increasing,  as  in  the  previous  test  cases,  an  unstable  limit  cycle 
would  exist. 

In  the  following  table  a  comparison  is  made  between  the  dynamic 
describing  function  curve  prediction  and  the  analog  computer  simula- 
tion of  the  amplitude  and  frequency  of  the  limit  cycle. 

Table  1.   Comparison  of  Limit  Cycle  Performance  Between  Dynamic 
Describing  Function  Curve  and  Analog  Computer. 


Saturation  Level 

Theoretical 

Computer 

For  N 

Amplitude 

uJ 

Amplitude 

uO 

+  1  volt 

1.6  volts 

1.07 

1.0  volt 

1.04 

+  2  volts 

2.9  volts 

1.1 

2.5  volts 

1.12 

+  3  volts 

2.7  volts 

1.27 

3.0  volts 

1.25 

+  5  volts 

3.2  volts 

1.28 

._ . 

3.0  volts 

1.26 

As  previously  noted  the  amplitude  of  the  initial  condition  varies 
in  the  vicinity  of  high  values  of  zeta^i.e.,  3       =  .6  to  .9.   The 
reason  for  this  unique  behavior  is  not  explainable  at  this  time. 
Possibly,  the  technique  used  in  this  thesis  has  reasonable  validity  for 
conditions  under  which  the  values  of  zeta  do  not  exceed  some  level,  say 
JT  -   .7 .      However,  is  all  cases  of  this  example  the  amplitude  of  A 
decreases  as  the  stable  limit  cycle  is  approached. 

Keeping  the  previous  limitation  in  mind  the  transient  response  was 
graphically  determined  for  an  initial  condition  of  5  volts.   This  value 
was  chosen  because  the  four  curves  behaved  the  same  for  this  initial 
condition  amplitude^!. e. ,  the  amplitude  of  A  decreases  as  the  limit 
cycle  was  approached. 
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In  Figures  4-24  to  4-27  are  the  calculated  transient  responses  and 
in  the  upper  right  hand  corner  of  each  plot  appears  a  brush  recording 
trace  of  the  actual  computer  response.   It  should  be  observed  that  the 
analog  computer  simulation  also  exhibits  a  stable  limit  cycle  for  the 
four  cases. 

4-3  Relationship  between  the  dynamic  describing  function  curve  and  the 

root  locus 

The  system  performance  of  a  nonlinear  system  depends  upon  the  motion 
of  the  characteristic  equation  roots.   In  using  a  graphical  technique 
which  assumes  second  order  dominance,  real  root  effects  are  neglected. 
Additionally,  one  would  suspect  that  real  roots  could  not  be  neglected 
for  a  fourth  order  system  since  there  are  regions  along  the  dynamic 
describing  function  curves  in  the  previous  sections  from  which  tangents 
can  be  drawn  to  the  J, ,  ■  1  curve  --  indicating  the  presence  of  real 
roots.   To  determine  the  values  of  the  characteristic  roots  at  discrete 
intervals  along  the  dynamic  describing  function  curve,  a  dynamic  root 
locus  plot  was  made. 

4-3-1     Dynamic  root  locus  plot 

Each  point  on  the  dynamic  describing  function  curve  can  be  identi- 
fied as  an  operating  or  M(  c*_  ,  (3  )  point  at  one  instant  of  time.  The 
roots  of  the  characteristic  could  be  obtained  by  inspection  for  each 
operating  point  since  zeta,  omega,  and  sigma  are  known.  Another  method 
is  to  substitute  the  coordinates  of  the  M(  at  ,  Q,  )  points,  N-  and  N  , 
into  the  general  characteristic  equation,  which  is  repeated  below,  and 
let  the  digital  computer  solve  for  the  roots. 

s4  +  7s3  +  14s2  +  (8  +  20N  )s  +  20NX  =  0  (4-1) 
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This  latter  technique  was  used.  When  the  roots  are  plotted,  we  obtain 
a  dynamic  root  locus  plot  which  determines  the  system's  characteristic 
roots  for  each  M(  <x_  ,  Q    )  point  shown  on  the  dynamic  describing  func- 
tion locus.   Dynamic  root  locus  plots  were  made  for  each  of  the  previous- 
ly determined  dynamic  describing  function  curves  from  thrt.  unstable 
point  to  the  linear  M-point  and  the  intermediate  solution  points. 

4-3-2     Summary  of  results 

Figures  4-28  to  4-33  are  the  dynamic  root  locus  plots  for  the  in- 
dicated examples.   Because  of  the  proximity  of  the  plotted  points,  the 

identified  roots  are  labelled  by  A.,  B.,  C,  and  D.  (i  =  1-4)  which 

J      l   i '   i*      l 

corresponds  to  the  unstable  condition,  limit  cycle,  an  intermediate 
value,  and  the  linear  M-point  respectively.   The  arrows  on  the  plots 
denote  the  direction  of  decreasing  values  of  A. 

In  all  cases  the  calculated  frequency  of  the  limit  cycle  points 
B.  agreed  accurately  with  the  dynamic  describing  function  curve  predic- 
tion.  However,  as  the  dominant  pair  of  conjugate  roots  moved  into  the 
left-half  plane,  the  real  root  effect  becomes  evident.   From  these  root 
locus  plots  the  form  of  the  transient  solution  can  be  determined  at  any 
instant.  When  two  pairs  of  conjugate  complex  poles  exist,  as  in  the 
vicinity  of  points  B.,  the  form  of  the  output  transient  solution  for  an 
initial  condition  input  is 

Q      (t)  =  Ae-l^t   cos  (  CJ  t  +  0)  +  Be  -  ST  t  cos  (  Uj  t  +  0  ) 

(4-2) 
If  the  second  term  is  much  smaller  than  the  first,  we  have  essentially 
second  order  dominance.   When  two  real  roots  and  a  complex  pair  of  con- 
jugate poles  exist,  as  in  the  vicinity  of  points  C,  the  transient  solution 
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Figure  4-30.   Dynamic  Foot  Locus  Plot  for  Figure  4-22, 
Curve  A. 
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Figure  4-31.   Dynamic  Root  Locus  Plot  for  Figure  4-22,  g 
Curve  B. 
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Fijjjire  4-152.      Dynamic   Root   Locus   Plot   for  Figure   4-22, 
Curve   C. 


fts 


is  of  the  form 

0O(t)  =  Ae-  ^H  cosO-Ot  +  0)  +  Be-  r,    t+De-  *i  t       (4-3) 
In  this  case  second  order  dominance  would  not  be  justified  unless  both 
real  roots,  (T .    and  <T  ,  were  far  enough  out  that  their  effect  is 
negligible  during  the  transient  response.   In  Figures  4-28  and  4-29 
one  of  the  real  roots  approaches  the  origin  for  decreasing  values  of  A 
which  accounts  for  the  real  root  response  from  the  analog  computer.   For 
the  other  cases  in  Figures  4-30  to  4-33  the  real  root  effect  is  not  ap- 
parent, because  the  system  always  approaches  a  stable  limit  cycle. 
The  mathematical  transient  solution  could  be  obtained  from  the  dynamic 
root  locus  plot  or  the  dynamic  describing  function  curve  by  solving  a 
series  of  differential  equations  -  a  tedious  and  unsatisfactory  method. 

The  reason  for  the  good  correspondence  between  the  graphical  techni- 
que and  the  analog  computer  is  that  several  seconds  have  elapsed  before 
the  real  root  is  noticeable.   By  this  time  the  transient  response  has 
decayed  to  a  low  value  such  that  only  a  small,  slowly  decaying  steady 
state  output  exists. 
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CHAPTER  V 

TWO  NONLINEAR  ELEMENTS  RELATED  TO  A  COMMON  SIGNAL  BY  A 

NONLINEAR  DIFFERENTIAL  EQUATION 

The  third  case  for  relating  two  nonlinear  elements  to  a  common 
reference  signal  is  more  complicated  than  the  previous  two  cases.   A 
typical  example  is  the  single  feedback  loop  system  shown  in  Figure  5-1 
where  the  input  signal  to  H     depends  upon  the  characteristics  of  the  non- 
linear element,  N  ,  when  both  elements  are  referred  to  the  common  signal 
excitation,  A.   Both  G1 (s)  and  G-(s)  are  assumed  to  have  low  frequency 
filter  characteristics  so  that  the  requirements  for  expressing  N  and 
N   in  terms  of  describing  functions  are  met.   This  system  will  first  be 
discussed  from  a  classical  viewpoint  and  lastly  related  to  the  para- 
meter plane  by  a  specific  example. 

5-1  Classical  approach 

In  reference  [4]  Schiring  discusses  an  approximate  method  for  deter- 
mining the  system  stability  of  nonlinear  systems  with  no  theoretical 
limitations  on  the  number  or  type  of  nonlinear  elements  as  long  as  each 
element  is  separated  by  an  adequate  low  frequency  filter.   However, 
when  the  nonlinear  elements  number  more  than  two  and  become  functions 
of  both  signal  amplitude  and  signal  frequency,  the  amount  of  trial  and 
error  calculation  becomes  prohibitive  so  that  direct  simulation  on  a 
computer  becomes  the  only  feasible  approach  to  establishing  system 
stability  and  limit  cycles. 

If  the  two  nonlinear  elements  depicted  in  Figure  5-1  are  amplitude- 
sensitive  only,  the  Bode  phase  diagram  can  be  used  to  find  the  resonant 
frequency,  CO  n»   The  critical  magnitude  criterion  may  be  used  to 
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Figure  5-1.   -Block  Diagram  of  Typical  Example  for 
Oase  3. 


graphically  determine  the  amplitude  of  the  limit  cycle  i.e., 

|«i«Uyl      |VA)I      lG2<J-J0)l      ■     N^A^ICOuio)!)''       (5-1} 

Two  additional  curves  above  and  below  this  critical  amplitude,  if  it 
exists,  establish  the  relative  stability  of  the  system. 

When  one  or  both  of  the  nonlinear  elements  in  this  problem  have 
hysteretic  or  memory-type  characteristics,  Schiring's  method  becomes 
exceedingly  tedious  because  trial  and  error  calculations  are  required 
to  meet  the  phase  and  magnitude  relationships.   While  Shiring's  method 
is  time-consuming,  it  provides  accurate  results. 

5-2   Parameter  plane  approach 

The  system  analyzed  in  the  preceding  section  is  now  analyzed  us- 
ing the  parameter  plane  method  which  first  requires  writing  the  system 
characteristic  equation.   As  assumed  before,  N1  and  N   are  any  two  gen- 
eral nonlinear  elements.   The  system  characteristic  equation  for  Fig- 
ure 5-1  is 

1  +  G1(s)N1G2(s)N2  =  0  (5-2) 

The  existence  of  the  N.N  product  may  impose  some  practical  restric- 
tions on  the  parameter  plane  method  because  as  noted  by  Schiring  when 
the  product  of  two  describing  functions  or  a  composite  describing  func- 
tion is  used  both  the  amplitude  of  the  input  signal  and  the  character- 
istics of  the  individual  nonlinear  elements  may  limit  the  accuracy  of 
the  results.   However,  there  are  no  mathematical  restrictions  on  the 
N-N-  product. 

When  N  and  N  are  amplitude-sensitive  only^i.e.,  not  functions  of 
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phase,  one  could  let  the  parameter  plane  variables  take  on  the  values 
<x-  =  N  N  and  Q>  -   K  the  system  gain,  for  example.   In  this  case  the 
system  is  mathematically  reduced  to  one  nonlinear  element  or  a  composite 
describing  function  where  the  M  -  point  variations  occur  in  the  ©<  or 
N  N?  axis  direction. 

For  the  more  general  case  N  could  be  an  amplitude-sensitive  ele- 
ment and  N.  a  memory  type  or  hysteretic  element.   In  this  case  character- 
istic equation  (5-2)  becomes 

1  +  G1(s)N1G2(s)(Na  +  jNb)  =  0  (5-3) 

where  N-  =  N  +  jN  has  been  expressed  in  terms  of  its  real  and  imagin- 
ary components.   Expanding  equation  (5-3)  results  in  the  expression 


1  4-  N1NaG1(s)G2(s)  +  jN1NbG1(s)G2(s)  =  0  (5-4) 

Letting  <^c  =  N  N  and  Q   =  N-N  ,  the  general  coefficient  term,  a  , 
takes  on  the  form 

ak  =  bk*.  +  jbk  §  +  dk  (5-5) 

where  the  b,  terms  are  the  constant  coefficients  of  the  real  and  imagin- 
ary components  of  the  nonlinear  elements,  and  the  d  variables  are  the 

constant  terms  only.   Equation  (5-5)  has  been  substituted  into  the  gen- 

v\ 

eral  characteristic  equation  summation  i.e.  ,  f(s)  =  2-  a   sk  and  solved 

7  p>o  k 

v 
for  zS-    and  Q  using  Chebyshev  functions  by  Siljak.   The  following 


solutions  apply: 


oc    =  .  Vl  +  B2D2 


2  2 

;1     +B2 
1  l  (5-6) 


'2D1  "  V2 


o    2    ,    „    2 
Bl     +B2 
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where 

\   -       I    ("Dk  bk  UJ  nk  Uk_l(  /  ) 
B2  -   £  (-Dk  bk  UJ  nk  Uk  (  I     ) 


(5-7) 


Dl  = 


1  (-D   dk  cjO  nK  Ukl(  /  ) 
D„  =   £_    (-Dk  d.   Ujnk  U,   (J  ) 


2 

The  U,  (  I    )  terms  are  Chebyshev  functions  of  the  second  kind.   Equa- 
tions (5-6)  can  be  readily  programmed  on  the  digital  computer  with  the 

aid  of  equations  (5-7)  to  obtain  curves  of  constant  I     and  (jj  .   How- 

n 

ever,  the  curves  of  constant  sigma  and  constant  omega  derived  in  sec- 
tion 2-1-2  which  use  the  binomial  coefficient  expansion  will  require 
modification  because  the  implicit  assumption  of  real  coefficients  is  no 
longer  valid. 

The  i    »  0  curve  obtained  from  the  above  relationships  can  be  used 
to  determine  the  amplitude  and  frequency  of  a  limit  cycle  as  well  as 

relative  stability.   Plotting  a  grid  of  the  variations  of  N.N   and  N..N,  , 

la      lb 

as  functions  of  A  and  '~0  ,  atop  the  parameter  plane  curves  will  deter- 
mine the  amplitude  and  frequency  of  the  limit  cycle  at  the  common  radian 
frequency  intersection  on  the  J    =  0  curve. 

5-3  An  example 

A  specific  example,  which  had  been  previously  worked  by  Donohue  [l]j 
using  Schring's  method  and  simulated  on  the  analog  computer  will  now  be 
reworked  using  the  approach  outlined  in  the  previous  section.   Figure 
5-2  is  a  block  diagram  of  the  experimental  system  where  N.  and  N  are 
saturation  and  backlash  elements  respectively  and  K  is  a  variable  gain. 
The  system  characteristic  equation  is 
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3     2 
s  +  3s  +  2  s  +  40  K  N. (N  +  jn  )  =  0  (5-8) 

la     b 

where  N  has  been  expressed  in  terms  of  its  real  component,  N  ,  and 
2  a 

imaginary  component,  N,  .   The  °C  variable  is  selected  as  N.N  and  the 
b  la 

£  variable  as  N.N,  . 
1  b 

The  variations  of  N  N  versus  N  N  with  respect  to  signal  amplitude, 
A,  and  radian  frequency,  Oj  ,  are  determined  with  the  aid  of  describing 
function  theory.   Referring  to  Figure  5-2,  the  magnitude  of  the  input 
signal  to  N   is  E  =  A   G.(juJ  )    while  the  magnitude  of  the  input  sig- 
nal to  N  is  Z  =  A  I  G  (j  00   )   N   I  G  ( j  oJ  )   .   From  the  literature, 
the  describing  functions  for  the  saturation  element,  N  and  the  back- 
lash element,  N„,  are  repeated  below. 


Nx  -  -^T      (sin"1  R  +  R  ->J  1-R2)  (5_9) 

N  =  Na  +  jNu  =  —  [  -^P-  +  -§-  +  sin  2&)  -  -*-  cos2G 
2         b   fT  \  4      2     ^f  /     Tr 

where  R  ■  1  and  ©=  sin    —  ,  c  being  the  value  of  the  projection  of 

E 
the  backlash  discontinuity  on  the  x-axis.   These  two  equations  are  read- 
ily programmed  on  the  digital  computer  to  obtain  a  grid  of  '-0  -A  curves 
for  use  with  the  parameter  plane  curves.   When  this  net  of  curves  is 
super- imposed  on  the  parameter  plane,  the  common  radian  frequency  on  the 
S    *  0  curve  will  determine  the  frequency  and  amplitude  of  the  limit 
cycle. 

The  values  of  K  in  equation  (5-8)  were  considered  in  reference  fl], 
namely,  K  =.15  and  K  =  .25.   Figures  5-3  and  5-4  show  the  2=0   curve 
for  K  =  .15  and  K  =  .25  respectively.   Since  the  other  f    -curves  are 
not  used,  they  are  not  shown.   The  stable  region  is  enclosed  by  the 
shaded  boundary  on  the  X  =   0  curve.   Along  the  jf  =  0  curve  are  indicat- 
ed the  0~)   values.   The  lighter  lines  shown  in  both  plots  are  the 
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variations  of  N  N  and  N  N  with  respect  to  amplitude,  A,  and  radian 
frequency,  6,)  .   The  magnitude  and  frequency  of  the  limit  cycle  is  in- 
dicated adjacent  to  the  common  radian  frequency  intersection.   Since  the 
magnitude  of  the  limit  cycle  in  reference  [1]  was  determined  as  the  in- 
put signal  to  the  backlash  element,  N  ,  it  was  necessary  to  multiply  the 
magnitude  of  the  limit  cycles  indicated  in  Figures  5-3  and  5-4  by 

I  G..(j  uj  )   N   I  G  ( j  oJ  )   ,  where  uj  is  the  limit  cycle  frequency,  to 
obtain  the  limit  cycle  magnitude  for  reference  purposes.   A  comparison  of 
the  amplitude  and  frequency  of  the  limit  cycle  between  the  parameter 
plane  approach  and  the  results  from  reference  [1]  for  K  =  . 15  and  K  = 
.25  are  tabulated  below.   The  first  two  rows  labelled  "Schiring's 
Method"  and  "Analog  Computer"  were  obtained  from  [1],   The  results  in- 
dicated in  the  last  row  were  obtained  from  the  parameter  plane. 


K  =  .15 

Amplitude 

'j  j 

Schiring's  Method 

.7  volts 

.51 

Analog  Computer 

.8  volts 

.49 

Parameter  Plane 

.61  volts 

.52 

K  =  .25 

Amplitude 

tj$ 

Schiring's  Method 

.9  volts 

.69 

Analog  Computer 

1.0  volts 

.65 

Parameter  Plane 

.82  volts 

.68 

It  should  be  noted  that  excellent  correlation  results  among  the 
three  methods.   The  parameter  plane  method  also  predicts  a  stable  limit 
cycle  because  the  signal  amplitudes,  A,  are  decreasing  as  the  3=0 
curve  is  approached  from  the  stable  region.   This  conclusion  also  agrees 
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with  the  theoretical  and  experimental  results  from  [1]. 

As  already  indicated  in  section  5-1,  no  attempt  will  be  made  to 
obtain  the  form  of  the  transient  response  since  the  necessarysigma-omega 
curves  using  complex  coefficients  are  not  available.   However,  the  present 
example  indicates  for  systems  containing  two  dependent  nonlinear  elements 
that  the  amplitude  and  frequency  of  the  limit  cycle  as  well  as  the  re- 
lative stability  can  readily  be  predicted  with  minimal  effort. 
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CHAPTER  VI 
CONCLUSIONS 

It  has  been  demonstrated  that  systems  containing  two  nonlinear 
elements  belonging  to  three  general  classes  can  be  conveniently  handled 
on  the  parameter  plane.  When  the  parameter  plane  technique  is  combined 
with  the  "'dynamic  describing  function"  curve,  three  useful  pieces  of  in- 
formation in  varying  degrees  of  utility  can  be  obtained  for  predicting 
system  performance. 

1.  The  dynamic  describing  function  curve  in  the  parameter 
plane  lying  between  0  —  Z    —  1  determines  the  system's 
stable  region. 

2.  When  the  dynamic  describing  function  locus  crosses  the 
3-0  curve,  one  can  determine  the  amplitude  and  fre- 
quency of  the  limit  cycle  as  well  as  whether  a  stable  or 
unstable  limit  cycle  exists. 

3.  From  the  motion  of  the  M-point  along  the  dynamic  describ- 
ing function  locus  the  transient  performance  can  be  close- 
ly determined  when  the  following  restrictions  are  satis- 
fied: 

a.  Second  order  dominancy  conditions  must  be  met. 

b.  The  radian  frequency,  GO  ,  and  the  damping  ratio,  t,    , 
change  slowly  with  time. 

As  in  any  investigative  study,  the  analysis  and  synthesis  of  non- 
linear systems  on  the  parameter  plane  is  far  from  completion.   From  this 
preliminary  study,  several  interesting  problems  arose  which  are  deemed 
worthy  of  further  investigation. 
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A.  Closed  loop  frequency  response 

As  already  noted  the  dynamic  describing  function  curve  describes 
the  motion  of  the  operating  point  during  the  transient  period.   In  refer- 
ence [7]  the  closed  loop  frequency  response  is  obtained  for  a  single  non- 
linear element  with  the  aid  of  the  Nichol's  chart.   A  question  that  im- 
mediately arises  is  the  possibility  of  obtaining  the  closed  loop  frequency 
response  of  a  system  containing  two  nonlinear  elements  from  its  dynamic 
describing  function  curve.   The  main  difference  is  that  the  dynamic  des- 
cribing function  locus  considers  a  series  of  points  rather  than  a  single 
point  at  a  time. 

B.  Transient  response  prediction 

The  graphical  technique  used  in  this  thesis  requires  second  order 
dominance  as  a  necessary  condition  for  practical  application.   Since  the 
possibility  of  real  roots  increases  greatly  with  higher  order  systems, 
a  modification  to  the  present  technique  which  incorporates  this  real 
root  effect  would  greatly  improve  the  accuracy  of  predicting  nonlinear 
transient  performance. 

C.  Constant  sigma-omega  curves  with  complex  coefficients 
The  recursion  relationships  used  to  obtain  the  constant  sigma- 
omega  curves  made  the  implicit  assumption  that  the  coefficients  of  the 
characteristic  equation  were  real.   Since  the  coefficients  become  com- 
plex when  hysteretic  or  memory  type  element  are  introduced  into  a  non- 
linear system,  an  immediate  requirement  is  that  the  previously  determin- 
ed recursion  relationship  be  modified  or  rederived. 
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APPENDIX  A 
GRAPHICAL  METHOD  FOR  DETERMINING  THE  TRANSIENT  RESPONSE 

A  graphical  technique,  developed  by  Siljak  [6],  was  used  for  deter- 
mining the  transient  response  of  a  system  from  its  dynamic  describing 
function  curve.   This  technique,  which  is  accurate  for  second  order 
dominance  conditions,  considers  the  amplitude  of  the  initial  conditions 
along  the  dynamic  describing  function  curve  as  forming  the  envelope  of 
the  transient  response.   At  each  M-point  the  graphical  value  of  i   and  cO 
from  the  zeta-omega  curves  or  <f  and  (0  from  the  sigma-omega  curves  are 
read.   The  zero  signal  crossings  are  determined  from  the  phase  relation- 
ship 0  ■   \   cxJ  dt.   For  an  initial  condition  input  the  first  zero  cross- 
ing is  assumed  to  occur  in  fr  /2  radians  and  subsequent  zero  crossings 
in  multiples  of  ^X   radians.   The  following  example  assumes  an  initial 
condition,  A  =  3  volts,  for  the  describing  function  curve  appearing  in 
Figures  4-2  and  4-3. 

A-l  Transient  response  from  zeta-omega  curves 

The  system  transient  response  can  be  determined  for  the  dynamic 
describing  function  curve  on  the  zeta-omega  curves  using  the  following 
method.   Initially,  a  curve  of  A  vs  j>  ,  Figure  A-l,  and  A  vs  oJ^  , 
Figure  A-2  are  plotted  using  the  amplitude  values  indicated  alongside 
the  dynamic  describing  function  curve  in  Figure  4-2,  and  the  values  of 
£      and  °\\    passing  through  these  points.   At  each  point  the  values  of 


uJyn   are  obtained  from  the  relationship  oJ  =  uJa  \|  /-  J 

From  these  curves  values  of  A,  t>   ,  and  UJ^  are  tabulated  in  Table 
2  incrementally  from  3  volts. 
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Table  2.   Data  for  Calculating  Transient  Responses  from  Zeta- 
Omega  Curves. 


A 

1 

C0_A 

hx 

JV 

J^ 

t 

AO 

3.0 

.3 

1.22 

.366 

2.73 

1.07 

0 

Al 

2.5 

.4 

1.29 

.516 

1.94 

1.08 

.45 

A2 

2.0 

.5 

1.5 

.75 

1.33 

1.3 

.85 

A3 

1.5 

.52 

1.76 

.915 

1.09 

1.5 

1.15 

A4 

1.0 

.475 

2.01 

.955 

1.05 

1.77 

1.5 

A5 

.8 

.435 

2.11 

.918 

1.09 

.9 

1.65 

A6 

0.0 

2.36 

The  reciprocal  of  1  ^^  gives  the  incremental  exponential  time  decay, 
T-  .   Starting  with  AO  =  3  volts,  the  first  exponential  decay,  approxi- 
mated by  ^  =  2.73,  is  assumed  to  extend  along  a  short  portion  of  the 
straight  line  until  Al  ~   2.5  volts.   Here  the  procedure  is  repeated  for 
point  Al  using  a  new  ^~  -  1.94.  When  the  response  reaches  zero  or  a 
steady  state  value,  the  plot  is  terminated  as  in  Figure  A-3.  At  each 
point  the  time  is  read  from  the  abscissa  and  tabulated  in  Table  A-l. 
Converting  '^-V  to  oJ  ,  a  plot  of  uJvs.  time  is  made  as  in  Figure  A-4. 
By  graphical  integration  the  zero  crossings  at  — r—  ,  Vt   ,  2  TV ,  etc. 
are  determined. 

The  final  plot  uses  the  amplitude  envelope  determined  from  Figure 
4-3  and  the  zero  crossings  from  Figure  A-4,  to  sketch  the  transient 
response  shown  in  Figure  A-5. 
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A-2  Transient  response  from  sigma-omega  curves 

When  sigma-omega  curves  are  available,  the  transient  response  can 
be  determined  more  rapidly  because  the  conversion  of  CO  to  u3*  and  vice 
versa  are  avoided  by  the  simple  expedient  of  using  constant  sigma  curves 
where  £^-&'   vice  constant  zeta  curves.   The  previous  example  is  reworked 
using  the  corresponding  sigma-omega  curves  appearing  in  Figure  4-3.   The 
indicated  values  of  A,  and  the  values  of  sigma  and  omega  passing  through 
these  points  are  used  to  plot  curves  of  A  vs  <T  and  A  vs  u3  which  appear 
in  Figures  A-6  and  A-7  respectively.   These  curves  are  used  to  tabulate 
the  following  table  for  values  of  A  selected  in  one-half  volt  increments 
from  three  volts. 

Table  3.   Data  for  Calculating  Transient  Response  from  Sigma- 
Omega  Curves. 


A 

x 

T- 

u) 

t 

AO 

3.0 

.35 

2.86 

1.18 

.0 

Al 

2.5 

.54 

1.85 

1.2 

.475 

A2 

2.0 

.75 

1.33 

1.3 

.85 

A3 

1.5 

.92 

1.09 

1.5 

1.175 

A4 

1.0 

.95 

1.05 

1.8 

1.52 

A5 

.5 

.9 

1.11 

1.98 

2.05 

A6 

.0 

3.16 

The  values  of  sigma  and  omega  are  read  from  their  respective  charts. 
The  exponential  decay  is  now  the  reciprocal  of  sigma  and  is  denoted  by 
the  column  marked  *>-  .   Commencing  with  A  =  3  volts,  the  envelope  is 
assumed  to  exist  along  the  tangential  line  terminated  by  'N-  ■  2.86 
until  A  -  2.5  volts  is  reached.   Time  is  read  along  the  x-axis  and  the 
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procedure  repeated  for  the  remaining  points  as  shown  in  Figure  A-8. 
To  obtain  the  zero  crossings  a  plot  of  u)  vs  t  is  made  (Figure  A-9) , 

and  by  graphical  integration  the  zero  crossings  are  determined.   The 
overall  transient  response  appearing  in  Figure  A-10  is  a  composite  plot 
of  the  signal  envelope  determined  in  Figure  A-8  and  the  zero  crossings 
in  Figure  A-9. 

A  visual  comparison  of  the  two  transient  responses,  namely,  Fig- 
ures A-5  and  A-10  shows  that  the  calculated  transient  responses  agree 
closely.   There  is  a  small  difference  in  the  terminal  end  point  of  the 
plots  but  this  is  an  insignificant  matter.   The  important  point  is  that 
both  methods  can  predict  the  form  of  the  transient  response  -  but  ac- 
curately only  under  second  order  dominancy  conditions.   In  addition, 
the  use  of  sigma-omega  curves  increases  the  rapidity  of  this  graphical 
technique  because  needless  conversions  are  eliminated. 
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APPENDIX  B 

DIGITAL  COMPUTER  PROGRAM  FOR  PLOTTING  CONSTANT 

ZETA-SIGMA-OMEGA  CURVES 

A  digital  computer  program  was  written  in  FORTRAN  60  for  use  on 

the  CDC  1604  computer  to  plot  the  curves  of  constant  zeta,  sigma,  or 

omega  used  in  this  thesis.   The  comment  cards  at  the  beginning  of  the 
program  indicate  how  to  use  the  computer  program. 
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